Week 02 - Aarhus Environs: Introduction to Reading and Plotting Spatial Data
Introduction into GIS in R
There are a few key spatial packages available for Spatial Analysis in R, which you need to install as you progress through these exercises. The most basic are
sffor working with vector datarasterandterrafor working with raster data
Task 1: Reading vector data
The sf package, created by Edzer Pebesma and colleagues,
has dramatically simplified reading vector spatial data into R.
In this exercise you will read in three shapefiles (parks,
playgrounds, and forests one point file and two polygon files) and one
geojson (shelters) using st_read(). If you’ve read in the
files correctly, you will see a standard R data frame except it will
show some header metadata about the file and you’ll see a special
geometry column which we will discuss later.
Instructions
- Load the
sfpackage. - All your datasets reside in the ‘data’ folder.
- Import the
forestsshapefile (“forests.shp”). - Import the
playgroundsshapefile (“playgrounds4326.shp”). - Import the
parksshapefile (“parks.shp”). - Import the
sheltersgeojson (“shelters.json”). - Use the
head()function and identify the first few features of each layer.
# Load the sf package
___(___)
# Read in the forests shapefile
forests <- ___("../data/forests.shp")
# Read in the parks shapefile
parks <- ___(___)
# Read in the playgrounds shapefile
playgrounds <- read_sf("../data/playgrounds4326.shp", crs = 4326)
# Read in the shelters json
shelters <- ___(___)
# View the first few features of all layers
___(___)Questions:
- How many features does each layer contain and what kind of geometry are they?
- What is the CRS value in these objects?
Solution
Reading layer `parks' from data source
`C:\Users\Adela\Documents\RStudio\1_Teaching\cds-spatial\data\parks.shp'
using driver `ESRI Shapefile'
Simple feature collection with 47 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 561480.2 ymin: 6210157 xmax: 581216 ymax: 6237475
Projected CRS: ETRS89 / UTM zone 32N
Reading layer `forests' from data source
`C:\Users\Adela\Documents\RStudio\1_Teaching\cds-spatial\data\forests.shp'
using driver `ESRI Shapefile'
Simple feature collection with 32 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 566833 ymin: 6212338 xmax: 582568.3 ymax: 6236083
Projected CRS: ETRS89 / UTM zone 32N
Reading layer `playgrounds4326' from data source
`C:\Users\Adela\Documents\RStudio\1_Teaching\cds-spatial\data\playgrounds4326.shp'
using driver `ESRI Shapefile'
Simple feature collection with 51 features and 5 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 9.990331 ymin: 56.03292 xmax: 10.30607 ymax: 56.26783
CRS: NA
Reading layer `shelters' from data source
`C:\Users\Adela\Documents\RStudio\1_Teaching\cds-spatial\data\shelters.json'
using driver `GeoJSON'
Simple feature collection with 19 features and 6 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 563459.7 ymin: 6212589 xmax: 578934.2 ymax: 6232898
Projected CRS: ETRS89 / UTM zone 32N
Simple feature collection with 6 features and 6 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 570961.1 ymin: 6216613 xmax: 577431.9 ymax: 6232898
Projected CRS: ETRS89 / UTM zone 32N
bookbar type navn skov
1 Ja Shelters Shelter i Gjellerup Skov <NA>
2 Ja Shelters Shelter i H\xf8rhaven i skoven <NA>
3 Ja Shelters Shelter i Lisbjerg gammel skov <NA>
4 Ja Shelters Shelter ved Moesg\xe5rd Strand <NA>
5 Ja Shelters Shelter i Mollerup Skov <NA>
6 Ja Shelters Shelter i H\xf8rhaven p\xe5 bakken <NA>
mi_style mi_prinx
1 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 1
2 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 2
3 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 3
4 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 4
5 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 5
6 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 6
geometry
1 POINT (570961.1 6223407)
2 POINT (576307.6 6219521)
3 POINT (572490.7 6232898)
4 POINT (577431.9 6216613)
5 POINT (574609.9 6229691)
6 POINT (576455.1 6219324)
Simple feature collection with 6 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 566833 ymin: 6223092 xmax: 576573 ymax: 6232903
Projected CRS: ETRS89 / UTM zone 32N
bookbar type navn skovnr areal_m2 min_x min_y max_x
1 Ja Skove_store Åbyhøj Skov 290 32656 -224407 191537 -224289
2 Ja Skove_store Brabrand Skov 240 39231 -225585 191149 -225221
3 Ja Skove_store Årslev Skov 250 76629 -228610 191124 -228087
4 Ja Skove_store Gjellerup Skov 190 216967 -224931 190847 -224289
max_y oprettet_a oprettet_d rettet_af rettet_dat
1 192017 az25000 2019-05-07 12:02:50.657 az25000 2019-05-07 12:02:50.657
2 191359 az25000 2019-05-07 12:02:50.657 az25000 2019-05-07 12:02:50.657
3 191374 az25000 2019-05-07 12:02:50.657 az25000 2019-05-07 12:02:50.657
4 191485 az25000 2019-05-07 12:02:50.657 az25000 2019-05-07 12:02:50.657
mi_style mi_prinx
1 Pen (2, 1, 32768) Brush (2, 32768, 16777215) 1
2 Pen (2, 1, 32768) Brush (2, 32768, 16777215) 2
3 Pen (2, 1, 32768) Brush (2, 32768, 16777215) 3
4 Pen (2, 1, 32768) Brush (2, 32768, 16777215) 4
geometry
1 MULTIPOLYGON (((571023.4 62...
2 MULTIPOLYGON (((570135.1 62...
3 MULTIPOLYGON (((567335.9 62...
4 MULTIPOLYGON (((571154.8 62...
[ reached 'max' / getOption("max.print") -- omitted 2 rows ]
Simple feature collection with 6 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 567350.2 ymin: 6210157 xmax: 581216 ymax: 6237475
Projected CRS: ETRS89 / UTM zone 32N
bookbar type navn areal_m2 oprettet_a oprettet_d
1 Ja Parker Aarslev Møllepark 40022 az25000 2019-05-07 12:02:19.407
2 Ja Parker Bananen 5667 az25000 2019-05-07 12:02:19.407
3 Ja Parker Bavneparken 35611 az25000 2019-05-07 12:02:19.407
4 Ja Parker Bækvejsparken 19846 az25000 2019-05-07 12:02:19.407
5 Ja Parker Engdalgårdsparken 81476 az25000 2019-05-07 12:02:19.407
6 Ja Parker Forteparken 65751 az25000 2019-05-07 12:02:19.407
rettet_af rettet_dat
1 az25000 2019-05-07 12:02:19.407
2 az25000 2019-05-07 12:02:19.407
3 az25000 2019-05-07 12:02:19.407
4 az25000 2019-05-07 12:02:19.407
5 az25000 2019-05-07 12:02:19.407
6 az25000 2019-05-07 12:02:19.407
mi_style mi_prinx
1 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215) 1
2 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215) 2
3 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215) 3
4 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215) 4
5 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215) 5
6 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215) 6
geometry
1 MULTIPOLYGON (((567616.2 62...
2 MULTIPOLYGON (((581047.6 62...
3 MULTIPOLYGON (((567728.9 62...
4 MULTIPOLYGON (((574082.9 62...
5 MULTIPOLYGON (((575090.4 62...
6 MULTIPOLYGON (((575844.2 62...
Simple feature collection with 19 features and 6 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 563459.7 ymin: 6212589 xmax: 578934.2 ymax: 6232898
Projected CRS: ETRS89 / UTM zone 32N
First 10 features:
bookbar type navn skov
1 Ja Shelters Shelter i Gjellerup Skov <NA>
2 Ja Shelters Shelter i H\xf8rhaven i skoven <NA>
3 Ja Shelters Shelter i Lisbjerg gammel skov <NA>
4 Ja Shelters Shelter ved Moesg\xe5rd Strand <NA>
5 Ja Shelters Shelter i Mollerup Skov <NA>
6 Ja Shelters Shelter i H\xf8rhaven p\xe5 bakken <NA>
7 Ja Shelters Shelter p\xe5 Vestereng <NA>
8 Ja Shelters Shelter i Moesg\xe5rd Skov <NA>
9 Ja Shelters Dagshelter i H\xf8rhaven <NA>
10 Ja Shelters Shelter i Brendstrup Skov <NA>
mi_style mi_prinx
1 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 1
2 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 2
3 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 3
4 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 4
5 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 5
6 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 6
7 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 7
8 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 8
9 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 9
10 Symbol (34,8388608,14,"MapInfo Real Estate",1,0) 10
geometry
1 POINT (570961.1 6223407)
2 POINT (576307.6 6219521)
3 POINT (572490.7 6232898)
4 POINT (577431.9 6216613)
5 POINT (574609.9 6229691)
6 POINT (576455.1 6219324)
7 POINT (573438.3 6227501)
8 POINT (576374.9 6217600)
9 POINT (576498.6 6219401)
10 POINT (572113.1 6227183)
Simple feature collection with 51 features and 5 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 9.990331 ymin: 56.03292 xmax: 10.30607 ymax: 56.26783
CRS: NA
First 10 features:
ogr_fid nr lokatin adresse
1 1 1 A. H. Wingesvej A. H. Wingesvej 13, 8000 Aarhus C.
2 2 2 Agernvej Agernvej 73, 8330 Beder
3 3 3 Baunevej Baunevej 15, 8361 Hasselager
4 4 4 Bellevue Fortevej 109, 8240 Risskov
5 5 5 Beringsminde Tranebærvej 44, 8220 Brabrand
6 6 6 Botanisk Have v. Peter Holms vej 17, 8000 Aarhus C.
7 7 61 Botanisk Have v. Mønsgade 12, 8000 Aarhus C.
8 8 7 Byvangen Grundvigsvej 55, 8260 Viby J.
9 9 8 Bækvejsparken Bækvej 31, 8340 Malling
10 10 9 Egå Engsø Viengevej 7, 8250 Egå
X_tmstm geometry
1 2020-02-14 00:00:00.0 POINT (10.20414 56.17545)
2 2020-02-14 00:00:00.0 POINT (10.20171 56.06382)
3 2020-02-14 00:00:00.0 POINT (10.09647 56.10748)
4 2020-02-14 00:00:00.0 POINT (10.24744 56.19133)
5 2020-02-14 00:00:00.0 POINT (10.12097 56.16046)
6 2020-02-14 00:00:00.0 POINT (10.19203 56.16155)
7 2020-02-14 00:00:00.0 POINT (10.19378 56.15858)
8 2020-02-14 00:00:00.0 POINT (10.17067 56.13048)
9 2020-02-14 00:00:00.0 POINT (10.19102 56.03292)
10 2020-02-14 00:00:00.0 POINT (10.22041 56.21385)
name long_name write copy
PCIDSK PCIDSK PCIDSK Database File TRUE FALSE
netCDF netCDF Network Common Data Format TRUE TRUE
PDS4 PDS4 NASA Planetary Data System 4 TRUE TRUE
VICAR VICAR MIPL VICAR file TRUE TRUE
JP2OpenJPEG JP2OpenJPEG JPEG-2000 driver based on OpenJPEG library FALSE TRUE
PDF PDF Geospatial PDF TRUE TRUE
MBTiles MBTiles MBTiles TRUE TRUE
BAG BAG Bathymetry Attributed Grid TRUE TRUE
EEDA EEDA Earth Engine Data API FALSE FALSE
OGCAPI OGCAPI OGCAPI FALSE FALSE
is_raster is_vector vsi
PCIDSK TRUE TRUE TRUE
netCDF TRUE TRUE FALSE
PDS4 TRUE TRUE TRUE
VICAR TRUE TRUE TRUE
JP2OpenJPEG TRUE TRUE TRUE
PDF TRUE TRUE TRUE
MBTiles TRUE TRUE TRUE
BAG TRUE TRUE TRUE
EEDA FALSE TRUE FALSE
OGCAPI TRUE TRUE TRUE
[ reached 'max' / getOption("max.print") -- omitted 62 rows ]
[1] "C:/Users/Adela/Documents/RStudio/1_Teaching/cds-spatial/Week02"
Well done, now you should see how easy it can be to read in shapefiles and you got your first taste of what an sf object looks like.
Task 2: sf objects are data frames
Spatial objects in sf are just data frames with some
special properties. This means that packages like dplyr can
be used to manipulate sf objects. In this exercise, you
will use the dplyr functions select() to
select or drop variables, filter() to filter the data and
mutate() to add or alter columns.
Instructions
- Load the
dplyrandsfpackages. - Use the
filter()function fromdplyron theparksobject to create a new data frame limited to parks greater than 5ha. - Use the
nrow()function on your new object to determine how many parks greater than 5ha are in the dataset. - Use the
select()function from dplyr to limit the variables in your over5ha dataset to justnavnandareal_m2and create a new data frame. - Use the
head()function to check which variables exist in your new data frame. Does the data frame only have thenavnandareal_m2columns (the answer is no, why)?
# Load the sf package
___
# ... and the dplyr package
___
# Use filter() to limit to over5ha parks
___
# Count the number of rows
___(over5ha)
# Limit to navn and areal_m2 variables
over5_lim <- over5ha %>% ___(navn, areal_m2)
# Use head() to look at the first few records
___(over5_lim)Solution
[1] 17
Great! You can see why the
sf package is so nice – your
spatial objects are data frames that you can smoothly manipulate with
dplyr. The number of parks over 5ha is 17 You may have
noticed that when you used select the default is to keep
the geometry column even if you didn’t explicitly list it as a column in
select.
Task 3: Geometry is stored in list-columns
A major innovation in sf is that spatial objects are
data frames. This is possible thanks, in part, to the list-column.
A list-column behaves, to a certain extent, like any other R column.
The main difference is that instead of a standard value such as a single
number, character or boolean value, each observation value in that
column is a piece of an R list and this list can be as complex as
needed. The list column allows you to store far more information in a
single variable and sf takes advantage of this by storing
all geographic information for each feature in the list.
In this exercise, you will convert the data frame to what’s called a
tibble with tibble::as_tibble() (Note that
dplyr::tbl_df() is now deprecated).
Instructions
- Load tidyverse in your workspace.
- Create a simple data frame
dfthat includes a single columnausingdata.frame(). - Add a list-column
bto your data frame with thelist()function. - Use
head()to look atdf. - Use
as_tibble()to convert the data frame to a tibble and print it to the console. This is just for cleaner printing. - Pull out the third observation from columns
aandbusingbaseR (you’ll need square brackets like[3]).
# Create a standard, non-spatial data frame with one column
df <- ___(a = 1:3)
# Add a list column to your data frame
df$b <- ___(1:4, 1:5, 1:10)
# Look at your data frame with head
___(df)
# Convert your data frame to a tibble and print on console
___(df)
# Pull out the third observation from both columns individually
df$___[___]
df$___[___]Solution
a b
1 1 1, 2, 3, 4
2 2 1, 2, 3, 4, 5
3 3 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
# A tibble: 3 × 2
a b
<int> <list>
1 1 <int [4]>
2 2 <int [5]>
3 3 <int [10]>
[1] 3
[[1]]
[1] 1 2 3 4 5 6 7 8 9 10
You now have a better sense of what a list column is. You can see how
it can be used to store far more information in a single variable than
other types of columns. These list-columns are how sf
stores detailed geographic information on each feature in a single
record. Converting the data frame to a tibble is not necessary but a
tibble can provide a better print out of the object.
Task 4: Extract geometric information from your vector layers
Sometimes you spatial data will come with attributes such as area, length, etc. Unless you have just created these yourself, you cannot be sure of their provenance and accuracy. They may have also changed in the course of your spatial transformations. If you want to know what is the actual area of your polygons, then you can extract it directly from the geometries.
There are several functions in sf that allow you to
access geometric information like area from your vector features. For
example, the functions st_area() and
st_length() return the area and length of your features,
respectively.
Note that the result of functions like st_area() and
st_length() will not be a traditional vector. Instead the
result has a class of units which means the vector result
is accompanied by metadata describing the object’s units. As a result,
code like this won’t quite work:
Instead you need to either remove the units with
unclass():
or you need to convert val’s class to
units, for example:
Instructions
- Check that your
sflibrary is still active and theforestsshapefile in memory - Compute the area of the forest units.
- Create a histogram of the areas using
hist()to quickly visualize the data spread. - Filter the forests object with
filter()and limit to forests withunclass(areas)> 200000 (areas greater than 20 hectares). - Can you plot the geometry and their names?
- Plot the geometry of the result with
plot()andst_geometry().
# Compute the areas of the forests
areas <- ___(forests)
# Create a quick histogram of the areas using hist
___(___)
# Filter to forests greater than 200000 (square meters)
big_forests <- ___ %>% ___(___(___) > 200000)
# Can you plot the big_forests with their names?(hint: check the plotting section in Task 5 below)
# Plot just the geometry of big_forests
___(___(big_forests))Solution
[1] "Gjellerup Skov" "Mollerup Skov"
[3] "Vestereng" "Tranbjerg Skov"
[5] "Vilhelmsborg Skov" "Skødstrup Skov"
[7] "Hørret Skov" "Brendstrup Skov"
[9] "Viby Høskov" "Fløjstrup Skov"
[11] "Kirkeskov" "Lisbjerg Skov, den gamle del"
[13] "Lisbjerg Skov, den nye del" "Thorsskov"
[15] "Hestehaven" "Skåde Skov"
[17] "Havreballe Skov" "Riis Skov"
[19] "Moesgård Skov" "Mariendal Skov"
Excellent! Computing geographic information for your vector layers
can be done with functions like st_area() and
st_length(). As you saw in this exercise, these functions
produce a result that can be used in additional calculations but you
need to be careful because the result is a units object
that requires a little additional processing like using
unclass().
Task 5: Plot vector spatial objects
The function for making a quick map/plot is a function you are
already familiar with, plot(). You can, for example, type
plot(my_data) to see your spatial object. The default,
though, may not be what you want. The plot() function, when
applied to sf objects, will create a set of maps, one for
each attribute in your data. Instead, if you want to create a map of a
single attribute you can extract that attribute using, as an example,
plot(my_data["my_variable"]).
Frequently you just want to plot the raw geometry with no attribute
color-coding (e.g., adding county boundaries to a map of points). For
this, you can use the st_geometry() function to extract the
geometry and plot the result. You can either create a new object or you
can nest st_geometry() within the plot()
function.
Often, you also want to plot multiple spatial objects together to see
if and how they relate. To do this, use the plot()function
twice in sequence, separated by ; with the second instance containing
the add=TRUE argument plot( ,add=TRUE). You
can use the col argument to differentiate each layer by
color.
Instructions
- Use
plot()to plot theforestdata using all defaults. - Plot just the
areal_m2attribute of forests. - Create a new object that is just the geometry of the forests object
with
st_geometry(). - Plot the geometry of the forests (the object you just created).
- Plot the geometry of the forests and the parks together, making the parks pink and the forests green.
Solution
Well done! Yes, these plots are not super pretty but you can’t beat
plot() for a quick look using few keystrokes. And remember
you can use plot(st_geometry(geo_object)) to plot just the
geometry of your object.
Task 6: Reading in raster data
The term “raster” refers to gridded data that can include satellite
imagery, aerial photographs (like ortophotos) and other types. In R,
raster data can be handled using the raster package created
by Robert J. Hijmans or by the newly developed terra
package by the same author. See which one you can get to work in your
workspace.
When working with raster data, one of the most important things to keep in mind is that the raw data can be what is known as “single-band” or “multi-band” and these are handled a little differently in R. Single-band rasters are the simplest, these have a single layer of raster values – a classic example would be an elevation raster where each cell value represents the elevation at that location.
Multi-band rasters will have more than one layer. An example is a color aerial photo in which there would be one band each representing red, green or blue light (RGB).
Instructions
- Load the
rasterorterrapackages. - Load the elevation grid
DNK_msk_alt.grdwith therasterfunction and assign toelevationobject - Read in the orthophoto image with
terra:rast()(this is multi-band raster called “Aarhus_1m.TIF”). - Use the
class()function to determine the class of each raster object you read in. - Use the
nlayers()andnlyr()function to determine how many bands/layers are in each object. - Use the
res()function to learn the resolution of the image
# Load the raster package
___(___)
# Read in the mound elevation single-band raster
elevation <- ____(___)
# Read in the orthophoto image multi-band raster
aarhus <- brick(_____)
# Get the class for the new objects
class(___)
class(___)
# Identify how many layers each object has
nlayers(___)
nlayers(___)
# Identify the resolution of each raster
___(elevation)
___(aarhus)Questions:
- What are the dimensions (number of columns and rows) of the
elevationraster? - What is the resolution of
aarhusorthophoto? How many layers does it contain and what do they represent?
Now you’ve learned how to read in single and multi-band rasters. You
should have noticed, based on the nlayers() function, that
the elevation object has a single layer and the
aarhus object has five (well, four useful ones, if you plot
it).
Solution
[1] "RasterLayer"
attr(,"package")
[1] "raster"
[1] "SpatRaster"
attr(,"package")
[1] "terra"
[1] 1
[1] 5
[1] 0.008333333 0.008333333
[1] 1 1
Task 7: Learn about your raster objects
Instead of storing raster objects in data frames, the
raster\terra package stores spatial data in specially
designed R classes that contain slots where the data and metadata are
stored. The data and metadata can be accessed using a suite of
functions. For example, the spatial extent (the bounding box) of the
object can be accessed with extent()\ext(), the coordinate
reference system can be accessed with crs() and the number
of grid cells can be determined with ncell().
Instructions
- You should have
rasterandterrapackages loaded, and have theelevationlayer and ortophoto image layer for Aarhus in memory. - Use the
extent()andext()functions to get the extent of theelevationandaarhuslayer. - Use the
crs()function to get the coordinate reference system ofaarhusandelevation. - Use the
ncell()function to determine how many grid cells are in theelevationlayer and theaarhuslayer.
# Get the extent of the elevation and aarhus object
___(___)
# Get the CRS of the aarhus and elevation object
___(___)
# Determine the number of grid cells in both raster objects
___(aarhus)
___(elevation)Solution
class : Extent
xmin : 8
xmax : 15.3
ymin : 54.5
ymax : 57.8
SpatExtent : 573334, 576668, 6223334, 6226668 (xmin, xmax, ymin, ymax)
[1] "PROJCRS[\"ETRS89 / UTM zone 32N\",\n BASEGEOGCRS[\"ETRS89\",\n DATUM[\"IRENET95\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]],\n CONVERSION[\"Transverse Mercator\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",9,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"easting\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"northing\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]]"
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +ellps=WGS84 +no_defs
WKT2 2019 representation:
GEOGCRS["unknown",
DATUM["Unknown based on WGS 84 ellipsoid",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1],
ID["EPSG",7030]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]],
CS[ellipsoidal,2],
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["latitude",north,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]
[1] 11115556
[1] 346896
Great work! Although rasters are not stored as data frames, the
metadata can easily be extracted using functions like
extent() or ext(), crs() and
ncell().
Question 3: It makes sense that the extents and number of cells for Danish national elevation data and Aarhus city would be different because their sizes diverge, but why are the units so different? What are these units?
Task 8: Plot your raster object
Similar to what you saw in the exercises related to vector objects
it’s often useful to quickly look at a map of your raster objects with
the plot() function.
The raster package has added useful methods for plotting both single
and multi-band rasters. For single-band rasters or for a map of each
layer in a multi-band raster you can simply use plot(). If
you have a multi-band raster with layers for red, green and blue light
you can use the plotRGB() function to plot the raster
layers together as a single image.
Instructions
- Plot the
elevationraster with theplot()function, it is a single-band raster. - Plot the
aarhusobject with theplot()function, it is a multi-band raster. - Plot the
aarhusraster withplotRGB()to see all layers plotted together as a single image. Check out the optionalstretchargument in this function, and set it tolinearto brighten the image up a little - Try to plot the two images together (on top of one another). Can you do it?
# Plot the elevation raster (single raster)
___
# Plot the aarhus raster (as a single image for each layer)
___
# Plot the aarhus raster as an image
____
# Plot the two images together (on top of one another). Hint: it only works if their CRSes are compatible. Return to this after Task 9 and 10 if the CRSes are different.
____Solution
Nice work! As you can see, the plot() function can be
used to plot single layers while the plotRGB() function can
be used to combine layers into a single image. Plotting two raster
objects of different extents, resolution and CRS can not be done without
additional wrangling
Task 9: Vector and raster coordinate systems
In order to perform any spatial analysis with more than one layer,
your layers should share the same coordinate reference system (CRS) and
the first step is determining what coordinate reference system your data
has. To do this you can make use of the sf function
st_crs() and the raster\terra function
crs().
When the geographic data you read in with sf already has
a CRS defined both sf and raster\terra will
recognize and retain it. When the CRS is not defined you will need to
define it yourself using either the EPSG number or the proj4string.
Instructions
- Ensure the packages
sfandraster\terraand the objectsforestsandplaygroundsandaarhusandelevationare loaded in your workspace. - Use
st_crs()to identify if a CRS exists and what it is for theplaygroundsandforestsobjects. - Use the
st_crs()function to define/assign a CRS to theplaygroundsobject, utilising the EPSG number 4326 as that is the original projection (now lost). - Use
crs()to identify if a CRS exists and what it is for theaarhusandelevationobjects. - Do not use the
crs()function to define/assign a CRS to theelevationobject as it already has a CRS of its own and renaming it won’t change the properties of the file.
Solution
Coordinate Reference System:
User input: ETRS89 / UTM zone 32N
wkt:
PROJCRS["ETRS89 / UTM zone 32N",
BASEGEOGCRS["ETRS89",
ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
MEMBER["European Terrestrial Reference Frame 1989"],
MEMBER["European Terrestrial Reference Frame 1990"],
MEMBER["European Terrestrial Reference Frame 1991"],
MEMBER["European Terrestrial Reference Frame 1992"],
MEMBER["European Terrestrial Reference Frame 1993"],
MEMBER["European Terrestrial Reference Frame 1994"],
MEMBER["European Terrestrial Reference Frame 1996"],
MEMBER["European Terrestrial Reference Frame 1997"],
MEMBER["European Terrestrial Reference Frame 2000"],
MEMBER["European Terrestrial Reference Frame 2005"],
MEMBER["European Terrestrial Reference Frame 2014"],
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[0.1]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4258]],
CONVERSION["UTM zone 32N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",9,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["Europe between 6°E and 12°E: Austria; Belgium; Denmark - onshore and offshore; Germany - onshore and offshore; Norway including - onshore and offshore; Spain - offshore."],
BBOX[38.76,6,84.33,12.01]],
ID["EPSG",25832]]
Coordinate Reference System: NA
# Assign the CRS to playgrounds - its CRS is actually 4326, it just got lost in
# transport
st_crs(playgrounds) <- 4326
playgrounds <- playgrounds %>%
st_set_crs(4326)
# Determine the CRS for the aarhus and elevation rasters
crs(aarhus)[1] "PROJCRS[\"ETRS89 / UTM zone 32N\",\n BASEGEOGCRS[\"ETRS89\",\n DATUM[\"IRENET95\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]],\n CONVERSION[\"Transverse Mercator\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",9,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"easting\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"northing\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]]"
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +ellps=WGS84 +no_defs
WKT2 2019 representation:
GEOGCRS["unknown",
DATUM["Unknown based on WGS 84 ellipsoid",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1],
ID["EPSG",7030]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]],
CS[ellipsoidal,2],
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["latitude",north,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]
# Assign the CRS to elevation
crs_2 <- "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"
crs(elevation) <- 4326Nice! You can determine what the CRS is using st_crs()
for vectors or crs() for rasters. If one doesn’t exist, you
can also use those functions to define the CRS.
Task 10: Transform your layers to a common CRS
In the previous exercise, when you ran st_crs() and
crs() you may have noticed that the CRS’ were different for
the different layers. The raster layer’s CRS began with
+proj=longlat and the vector layer’s began with
+proj=utm. In order to use these layers together in spatial
analysis we will want them to have the same CRS.
In this exercise you will transform (sometimes this is called
“project”) the objects so they share a single CRS. It is generally best
to perform spatial analysis with layers that have a projected CRS (and
some functions require this). To determine if your object has a
projected CRS you can look at the first part of the result from
st_crs() or crs() – if it begins with
+proj=longlat then your CRS is unprojected.
Note that you will use method = "ngb" in your call to
projectRaster() to prevent distortion in the elevation
image.
Instructions
- Use the
crs()function with argumentasText = TRUEon theaarhuslayer to get the CRS as a string and save this asthe_crs. - Use
st_transform()to transform the vectorplaygroundsobject to the CRS inthe_crs. - Use
projectRaster()to transform the rasterelevationobject to the CRS inthe_crs. This will take a few seconds. - Use
st_crs()onplaygroundsandelevationto confirm that they now have the same CRS. They should all have a CRS that starts with+proj=utm.*
Solution
# Get the CRS from the aarhus object
the_crs <- crs(aarhus)
# , asText = TRUE)
# Project playgrounds to match the CRS of aarhus
playgrounds_crs <- st_transform(playgrounds, crs = the_crs)
# Project to match the CRS of aarhus
elevation_crs <- projectRaster(elevation, crs = the_crs, method = "ngb")
# Look at the CRS to see if they match
st_crs(playgrounds_crs)Coordinate Reference System:
User input: PROJCRS["ETRS89 / UTM zone 32N",
BASEGEOGCRS["ETRS89",
DATUM["IRENET95",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]],
CONVERSION["Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",9,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
wkt:
PROJCRS["ETRS89 / UTM zone 32N",
BASEGEOGCRS["ETRS89",
DATUM["IRENET95",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]],
CONVERSION["Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",9,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
Coordinate Reference System:
Deprecated Proj.4 representation:
+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
WKT2 2019 representation:
PROJCRS["ETRS89 / UTM zone 32N",
BASEGEOGCRS["ETRS89",
DATUM["IRENET95",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]],
CONVERSION["Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",9,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
Great work! This may be the least fun part of spatial analysis but knowing how to do it will save you from a lot of frustration later.
Task 11: Plot vector and raster together
If the layers do not share a common CRS they may not align on a plot.
To illustrate, in this exercise, you will initially create a plot with
the plot() function and try to add two layers that do not
share the same CRS. You will then transform one layer’s CRS to match the
other and you will plot this with both the plot() function
and functions from the tmap package.
Note that for this exercise we returned all the layers to their original CRS and did not retain the changes you made in the last exercise.
With the plot() function you can plot multiple layers on
the same map by calling plot() multiple times. You’ll need
to add the argument add = TRUE to all calls to
plot() after the first one and you need to run the code for
all layers at once rather than line-by-line.
Instructions
- Try to plot the
playgroundsobject on top of theaarhusobject withplotRGB(aarhus)followed byplot(playgrounds, add = TRUE). Do you see the playgrounds? If not, use thest_transformto project to shared CRS. - Re-run the
plot()code from the instruction above, color theplaygroundsgreen and make the points stand out larger by setting thelwdargument to 3. You can also increase the transparency ofaarhus, by setting itsalphaat around 100. - Run the given
tmapcode. Note thattm_rgb()is used for multi-layered raster
Solution
# Plot aarhus and playgrounds (run both lines together) Do you see the
# playgrounds?
plotRGB(aarhus, stretch = "lin") #, alpha = 100)
plot(playgrounds_crs, col = "green", lwd = 3, add = TRUE)[1] "Modes: list, character"
[2] "names for target but not for current"
[3] "Attributes: < Modes: list, NULL >"
[4] "Attributes: < names for target but not for current >"
[5] "Attributes: < current is not list-like >"
[6] "Length mismatch: comparison on first 1 components"
# Run the tmap code
library(tmap)
tm_shape(aarhus) + tm_rgb() + tm_shape(parks) + tm_polygons(col = "green") + tm_shape(playgrounds_crs) +
tm_dots(col = "yellow", size = 1)
Great work! As you noticed, you mostly can’t plot layers together if
they don’t have the same CRS. You’ll see later that there are exceptions
but it is definitely best practice to ensure the layers you’ll work with
have a common, projected CRS.
Task 12: The Detective
I have prepared a mystery.zip folder for you at this
link: Download the mysterious
files. They are all in Denmark, but they have different CRSs and
file formats. Your job is to read them in R and create two plots: 1) one
where you make all layers line up and 2) where you select and filter
layers to create as pretty map as
If you need additional layers, check the Bonus section
below
Bonus: Loading data from the internet
There are a lot of online resources for spatial data, such as OSM,
AirBnB, GoogleMaps, DivaGIS, and GADM databases. Load some of this
spatial data - vector and raster - using the geodata
package. GADM has data for all current countries and is well-integrated
with R. You can load the data directly into R with the function
geodata::gadm() and start processing. The available data
are:
- SRTM 90 (elevation data with 30s - 3s resolution between latitude -60 and 60)
- World Climate Data (Tmin, Tmax, Precip, BioClim)
- Global adm. boundaries (different levels)
In the case of GADM you must also provide the level of administrative subdivision (0=country, 1=first level subdivision)
Instructions
- To read elevation with the
elevation_30s()function from thegeodatapackage, select the following arguments:countryneeds a 3 letter ISO code for Denmark; getData(‘ISO3’) let’s you see the codes, you can also specifyLatandLongarguments instead of country, if you need a specific area that overlaps multiple countries.maskneeds to be set to TRUE as we wish to mask surrounding countries.
- to read administrative data with
gadm()function, choose the following:countryneeds a 3 letter ISO code for Denmark; getData(‘ISO3’) let’s you see the codes,levelshould be set to the level of administrative subdivision (0=country border, 1=first level subdivision, eg. Midtjylland region, 2 = municipalities,… ).